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Abstract 

A partial fraction decomposition of the Fermi function resulting in a finite sum over 
simple poles is proposed. This allows for efficient calculations involving the Fermi 
function in various contexts of electronic structure or electron transport theories. The 
proposed decomposition converges in a well-defined region faster than exponential and 
is thus superior to the standard Matsubara expansion. 
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1. Introduction 

Many problems in electronic structure and electron transport calculations involve 
the evaluation of integrals containing the Fermi function. These are in general dif- 
ficult to compute and therefore several approximation schemes have been developed 
U 12. HI H S @] . Among them the Sommerfeld expansion |1] and the Matsubara ex- 
pansion |12|] being the most prominent ones. While the former is by construction useful 
for low temperatures, the latter provides in principle a way to cover the range from 
low to high temperatures. Moreover, it turns out that the expansion in a (finite) sum 
of simple poles is particularly suitable for evaluating the integrals by making use of 
the residue theorem. For example, finite temperature charge density calculations only 
require the evaluation of a Green's function at a finite set of energies (0, given by 
the poles of the expansion. Recently the same concept was used for the auxiliary den- 
sity matrix propagation in the context of time-resolved electron transport in molecular 
wires |@]. The major disadvantage of the Matsubara expansion consists in its poor con- 
vergence behavior, the error decreasing only linearly with the number of terms in the 
expansion. Here we derive an expansion of the Fermi function in terms of simple poles 
with particularly simple coefficients. We will show that it converges very rapidly with 
increasing order of the expansion in a well-defined region which is found to increase 
linearly with the order. 

For the following discussion it is convenient to write the Fermi function /(e) in 
terms of a dimensionless variable x, 

/w = t:^t x = (1) 

1 + e-^ kT 



Preprint submitted to Journal of Computational Physics 



March 27, 2009 



where yu is the chemical potential, T is the temperature, k is the Boltzmann factor and s 
denotes the energy. The expansion consists in finding a partial fraction decomposition 
with simple poles of the form 

2 ^—^ X - X,, 

where Ap are expansion coefficients and Xp are (possibly complex) poles. For practical 
purposes the sum over p is truncated and the Fermi function is approximated by f{x) ^ 
/iv(x) with being the number of terms in the expansion. 

For example, the well-known Matsubara expansion |2] is given in terms of the 
purely imaginary zeros x„ of the denominator in Eq. ([T]l, Xp = m{2p-l), which yields 
coefficients Ap = 1 and gives 



X ---V I ^ + 1 ] 



(3) 



For — > oo in Eq. (|3]l the expansion becomes exact. However, the convergence is very 
slow, which renders the application of this expansion impractical especially for low 
temperatures. 



2. Partial Fraction Decomposition 

The proposed partial fraction decomposition (PFD) is obtained by firstly writing 
Eq.Oas 1 5] 

1 1 1 sinh(x/2) 

m = T - T tanh(^/2) = - - , (4) 

2 2 2 2cosh(x/2) 

and secondly by expanding numerator and denominator in a power series, truncating 
the respective sums, such that the degree of the polynomial in the denominator is larger 
than the degree of the numerator polynomial. This procedure gives 

1 1 Pn-i(x/2) 

2 2 Qn(x/2) 
with polynomials 

^ ^2m+l N ^2m 

m={) m={) 

This construction allows for a PFD, i. e. an expansion of the form 

Pn-i(x/2) ^ A / Ap ^ Bp \ 
Qn(x/2) ^ \ x/2 - Xp x/2 + Xpj' 

Here, +Xp are the zeros of the polynomial Qn, which appear in pairs since Qn contains 
only even powers of x/2. It can be shown that the zeros can be obtained as Xp - -^z^, 
whereby the Zp are the eigenvalues of the following matrix Jsl, 

Zij = 2ii2i-l)6jj+i-2Ni2N-l)6iN, i,j^l,...,N. (8) 
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The eigenvalues can be efficiently calculated using standard methods. In Fig.[T]we have 
plotted the poles of Eq. (|7]i, given as +2xp= + 2 -y^, for three sets of eigenvalues Zp for 
different sizes A^xA^ of the matrix Z,y. As can be seen from the figure some of the poles 
+2xp arising from the PFD are purely imaginary and are close to the Matsubara poles. 
On the other hand there are also poles with a non-vanishing real part, which display 
an irregular distribution. These very poles improve considerably the approximation for 
the Fermi function as we show below. 




Refxj Re(x) Re(x) 

Figure 1: Poles (symbol •) of the PFD expansion, i.e. ±2-sJzj^, with Zp the eigenvalues of matrix (s) for 
vaiious orders N. For compaiison the purely imaginary poles (symbol x) of the Matsubara expansion (5) are 
shown as well. Note the different scales of the three graphs. 



It remains to determine the corresponding expansion coefficients Ap and Bp in 
Eq. (|7|l- Multiplying both sides of this equation by (x/2 - Xk) and letting x — > 2xk 
leaves on the right side of Eq. only the term Ak, which is thus given as 

Ak = hm (x/2-Xk) = lim ^ . ^ . ■ (9) 

x^ixi, Qn(x/2) n^o QNixk+rj) 

By means of the definitions ^ one finds from this limit Ak = \ and similarly Bk = I. 
Thus we arrive at the main result of this paper: The Fermi function can be approximated 
by the finite sum 

fN{x) = ^ - y (— ^ + — ^) . (10) 

2 j-[\x + 2^jrp x-2^jrpj 

with Zp the eigenvalues of matrix ([8]). The formal structure of this approximation is 
similar to the Matsubara expansion (|3j. However, taking advantage of having complex 
rather than purely imaginary poles makes the PFD for given order of the expansion 
vastly superior to the Matsubara expansion. This can be seen in Fig.|2] where we have 
shown both expansions for different orders A^. Whereas the Matsubara expansion (O in 
Fig.|2^ does not give a reasonable representation for any of the orders shown, the PFD 
expansion ( fTOl i in Fig.|2]3 improves rapidly with increasing order 



3. Convergence properties 

From Fig.|2] it becomes clear that the PFD is indeed converging faster than the 
Matsubara expansion. In the following we will quantify the rate of convergence as 
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Figure 2: Approximated Fermi function /ai(jc) for expansion orders N = 2,8, 32, 128 (full lines). Panel a: 
Matsubara expansion Eq. ^3), panel b: partial fraction decomposition Eq. jlO) . The curves are shown for 
X < only and are vertically shifted by 0.2 for better visibility. The exact Fermi function Eq. {T) is denoted 
by dotted lines. 



— > oo and give a range for x where this convergence behavior can be expectedQ To 
this end we define the deviation of the finite expansion from the exact function as 

6Mx) ^ fix) - Mx). (11) 

Regarding the PFD one makes two observations: First, in terms of the scaled variable 
y = x/4N one finds in the limit of large A^, 

f for > -1 

lim ^/.(x.,4A.) = I t (1 , 4^) ^ : (1 , i) , (12) 

The asymptotic function ( fT2b is shown along with deviations dfyiy) for various finite 
in Fig.[3h. Second, in the range -4A^ < x, i. e. -1 < y, the rate of convergence is 
given by the asymptotic expression 

SMx^ym - = (13) 

which due to the factorial in the denominator decreases faster than exponential. Eqs. ( fT2b 
and ( fT3] l are the main results of this section. They corroborate the statement that the 
PFD is expected to yield a better convergence and allow to estimate the error in actual 
calculations. 

In the remaining part of this section we will justify and discuss Eqs. iT2\ and ( fT3T l. 
Considering the case y < -I, one finds from Eq. ( flOt . that a finite expansion behaves 
as fN(x) oc 1 /x for x — > -oo. Since the Fermi function gives f(x) = 1 for x — » -oo one 
expects qualitatively the behavior given in Eq. (fT2l i. This holds true for any expansion 



In this section we consider only negative arguments x < 0, but the discussion applies in an analogous 
manner also to x > 0. 
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Figure 3: Deviation 6fN of tlie approximated Fermi function form the exact one as defined in Eq. 11 U . Left 
panel: as a function of the scaled argument y = xjAN for N = 2, 8, 32, 128 (solid lines) and the 
asymptotic expression (dashed line) as given by Eq. )12t . The dotted line shows r5/i28 for the Matsubara 
expansion. Right panel:i5/A> as a function of the expansion order N for three arguments x = -5, -25, -125. 
We compare the Matsubai'a expansion (dotted fines) and PFD (solid lines). For the latter case we show also 
the asymptotic behavior according to Eq. )13) by dashed lines. 



resulting in a finite sum over simple poles including the Matsubara expansion, which 
is shown for = 128 as dotted line in Fig.[3h. In order to verify that this behavior is 
indeed restricted to y < -1, or equivalently to x < -4N, we write the polynomial 
from Eq. (|6]l explicitly as 



N 



2m 



QN{xl2)^QN{y2N)^ } q,„N{y) with ^„,^(3;) = l-li-/-. (14) 



m=Q 



(2ot)! 



Assuming y < -1, we see that the ratio of two successive terms 

q,„N(y) 2 (^Nf 



qm-iNiy) 2m(2m-l) 

is always larger than 1 ; the terms are monotonically increasing. Thus terms with m » 1 
dominate the sum and we replace the coefficients in q„N by the coeflicient from q^^, 
i.e. instead of the sum (fl4l l we define 



v-i (2N)^-' o 

e^(y2A^) = 2]§"'^(3') with q„,^(y)^'-—L_y2>n^ (16) 

It turns out that in the limit N oo this sum becomes equal to Qf^{y2N), which can 
be seen by considering the difference of the newly defined terms in Eq. ( fT6l ) from the 
original terms in Eq. (fT4l i. For m - N - n one gets 



ciN-nAy) (2N-2n)\ N 

This expression vanishes for n - Q and can be made arbitrarily small by increasing 
for all n ^ N. Terms with larger n can be neglected because they are exponentially 
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small compared to those with smaller n. Since the sum ( fTSI l is a geometric series we 
obtain 

„2 



QN(y2N) ^ QM(y2N) ^ qMNiy) for N ^ (18) 



/-I 

Analogous considerations for the other polynomial from Eqs. ^ yield 



with pmNiy) - (3'2A^)^'"^V(2ot+1)! and we get as an approximation for the ratio, once 
again using » 1 , 

PN-i{y2N) 1 ^ 4A^ 

ew(3;2A^) ~ y x' ^ ' 

which explains the asymptotic behavior of 6 fN(x) in Eq. ( fT2b for x < -4N. 

Turning now to the case y > -1, we firstly note that there is a crossover for the 
ratio (flST l at \y\N; whereas for m < \y\N the terms are increasing, they decrease for 
m > \y\N. Thus, for large the polynomial expression (|5]l converges to the exact 
expression (|4|l and the deviation dfyix) vanishes as given by Eq. (fT2l) for x > -AN. In 
order to quantify the rate of convergence it is useful to define complementary sums to 
Pn and Qm, namely 

mts',, (2m+l)! J^^^ (2m)! 

Therewith the deviation reads 

sinh(3/2A^) - PN-\iy2N) sinh(y2A^) 



6MyAN) 



2 cosh(y2A?) - 22a'(3'2A?) 2 cosh(3;2A?) 
^ (e^(y2A?) - PN-iiy2N)) . (22) 

The approximation in the second line applies to large values of A^. We can choose, for 
any given y < Q, N sufficiently large such that exp(-3'2A^)» exp(+3'2A^). The infinite 
sums defined in (ISTT l become small compared to the exponentials, QN{y2N)^\^eK^{—y2N) 
and PAr_i(3;2A^)<Kl<K exp(-3'2A^). It remains to estimate their behavior for large 
which can be done in analogy to the considerations for and P^, cf. Eqs. (fTSl l and 
(fT9] l. Here the ratio of successive terms as defined in Eq. ( fTSI l is always smaller than 1 
and the first terms in the sum can be used to estimate the sums. One gets 

PN{y2N) ^ -^PNAy) and QN{y2N) ~ -^qN+iMy)- (23) 

This directly leads to Eq. ( fTsT l and concludes the derivation. 

Fig.[3j' shows this estimate along with the numerically calculated deviation dfy as 
a function of the expansion order A^ for selected values of x. Even for small values of 
A^ an overall good agreement is found. Moreover, one sees that the deviation Sfy is of 
order 1 as long as N < -x/4. However, for A^ > -x/4 (this is where the dashed lines 
start) it decreases very rapidly due to the factorial in the denominator in Eq. (fTST l. 
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4. Conclusions 



We have proposed the expansion (fTOl ) of the Fermi function ([TJ by using a partial 
fraction decomposition. Its appHcation requires only the diagonalization of a matrix, 
given in Eq. which has the same dimension as the expansion. The expansion 
converges faster than exponential with increasing order for arguments \x\ < AN. In 
other words, the approximation becomes not only more accurate for higher orders, it 
can also be used for a wider range of arguments. An estimate for the error is explicitly 
given by Eq. ( fT3l l. Due to the beneficial convergence properties and the straightforward 
implementation we expect the PFD to be of great value in any application based on an 
expansion of the Fermi function as sum of over simples poles. Finally, we would like 
to notice that an analogous expansion can be found for the Bose-Einstein distribution. 
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